Beam propagation management (BPM) is always a renewing topic and has attracted a great deal of interest, due to its potential applications in all-optical steering. As a possible route of BPM, spatial solitons have been systematically studied in different physical settings.^[1–4] However, in most of the literature, researchers shed their light on the transverse modulation of the employed potentials. Usually, by modulating optical mediums along both transverse and longitudinal directions, we can construct the so-called modulated photonic lattice, which can spawn many interesting phenomena that do not occur in purely transverse configurations, such as optical Bloch oscillations,^[5] Zener tunneling and Rabi oscillations,^[6,7] adiabatic light transfer,^[8] photon bouncing ball,^[9] resonant suppression of light-coupling,^[10] defect-free surface waves,^[11] diffraction-managed solitons,^[12,13] and light bullets,^[14] among many others (see Ref. [15] and references therein for a comprehensive review of modulated photonic lattices).

In the last decades, there have been widespread efforts for exploring the fractional effects, e.g., fractional Talbot effect,^[16] and fractional Josephson effect,^[17] which originate from the fractional Schrödinger equation (FSE). As a generation of the ordinary Schrödinger equation (SE), the FSE is developed by a straight substitution of fractional spatial derivatives for the usual ones.^[18] Although the concept of the FSE has received growing attention these days, most of the current research focuses on the mathematical aspects of theory and the evolution behavior of wave packets in simple potentials.^[19] Until 2015, by considering the formal equivalence between the SE and the optical wave equation, Longhi proposed a feasible scheme for realizing the fractional quantum harmonic oscillator optically.^[20] Following this pioneering work, a series of intriguing works based on the linear or nonlinear FSE have been reported. In the linear optics domain, zigzag trajectory propagation,^[21] PT-symmetry,^[22] diffraction-free beams,^[23] optical Bloch oscillation and Zener tunneling,^[24] beam propagation management in double-barrier potential,^[25] Anderson delocalization of light,^[26] and Airy trajectory engineering^[27] have been studied. Meanwhile, gap solitons,^[28] parity-time solitons,^[29] surface solitons,^[30] vortex solitons,^[31] and Hermite–Gaussian-like solitons^[32] modeled by the nonlinear FSE have also been investigated.

We notice that by introducing the longitudinal modulation of the transverse potential (Gaussian or periodic), resonant mode conversions and Rabi oscillation can be realized.^[33] When the fractional diffraction term is modulated longitudinally, one may observe periodic oscillation behavior of a Gaussian beam.^[34] Despite the above progresses, there still exists an open problem: What are influences of longitudinal modulations of the linear index potential on the propagation of a chirped Gaussian beam? which is the central goal of this work.

We start our analysis by investigating the propagation dynamics of a light beam modeled by the FSE as follows:

At the beginning, we consider the limiting case α = 1, which is solvable. If d(ξ) ≡ 0, equation (1) will be degenerate to the case addressed in Ref. [23]. When C = 0, equation (10) admits an approximated solution

3.1. Case 1: d(ξ)≡ 1

We firstly consider the simplest situation: C = 0. Under this condition, according to Eq. (9), the trajectory of all the frequency components is κ (k;ξ) = k – ξ/2, as shown in Fig. 1(a). In order to acquire the analytical solution to Eq. (10), we introduce an key parameter

which represents the accumulated phase. For those frequency components k ⩽ 0, k – ξ/2 < 0 always holds regardless of ξ, so the exact expression of the phase is ϕ(k;ξ) = –kξ + ξ²/4. For k > 0, k – ξ/2 is positive in the interval ξ ∈ (0,2k).

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (a) The trajectories of three different frequency components κ (k;ξ) = k – ξ/2. Solid blue: k = km (km is the maximum frequency component); dashed blue: k = –km; solid orange: k = 0. The horizontal line κ = 0 and line k = km intersect at the shift point (2km,0), after which κ(ξ) < 0 always holds. (b) Dependence of the integral form-factor on the linear chirp.

In the interval ξ ∈ (2k,∞), k – ξ/2 is negative, thus ϕ (k;ξ) satisfies the following equation:

As we all know, the beam’s profile is dictated to a large degree by the low frequency components. Therefore, for long-distance propagation, 2k² is far less than −kξ + ξ²/4, which means that we can make a reasonable approximation, i.e., ϕ(k;ξ)≈ –kξ + ξ²/4. Therefore, the analytical solution to Eq. (1) in real space can be written as

We can deduce the trajectory of the beam in real space as follows:

To verify the above analytical analysis, we have directly simulated the evolution behavior of the Gaussian beam by employing the so-called split step Fourier method. Due to the existence of linear index potential, the Gaussian beam depicted in Eq. (11) is mainly in the κ < 0 region [Fig. 2(b)], so that the beam will travel in a straight line [Fig. 2(a)], which is same as the C ≠ 0 case modeled by the FSE without a potential.^[23] All the observed phenomena prove that good agreements are obtained between the numerical results and the theoretical predictions. However, the integral form-factor, which is defined as and manifests the localization degree, is an exception. Concretely, according to its definition and Eq. (13), the integral form-factor should be a constant. Unfortunately, as illustrated in Fig. 2(c), it decreases first and then becomes a constant at ξ ≈ 3.32, which contributes to the splitting term kξ – ξ²/4.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Evolution behavior of a Gaussian beam without chirp in the FSE: (a) in (η,ξ)plane, (b) in (κ,ξ) plane. The dashed curves are analytical predictions. (c) The dependence between localization degree and propagation distance. In all the panels, α = 1, C = 0, η0 = 10, and σ = 0.25.

If C < 0 and | C| > k_m (k_m is the maximum frequency component), in both κ and η space, the beam will exhibit similar evolution behavior as the C = 0 situation [Figs. 3(a) and 3(b)]. However, it is worth mentioning that the integral form-factor under this condition is independent of ξ, as shown in Fig. 3(c). This is because k + C – ξ/2 is always negative, and then Eq. (11) can be calculated exactly as ϕ(k;ξ) = –kξ – Cξ + ξ²/4.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The same as Fig. 2, but with (a)–(c) C = –15 and (d)–(f) C = 15. The horizontal line in (e) corresponds to κ = 0. Those dotted vertical lines in (d)–(f) annotate the position of shift points at which κ changes its sign.

For C > 0, by taking the central frequency as an example, one may find its sign changes at ξ = 2C, thus following the same procedure as above, we can obtain the approximate expression for Eq. (11) as follows:

and the corresponding trajectory in real space is

which is consistent with the numerical simulation of the propagation displayed in Fig. 3(d). Meanwhile, in the κ space, the theoretical trajectory agrees well with the numerical simulation [Fig. 3(e)]. However, through an overall consideration of the effect of linear chirp on the dynamics characteristics of a Gaussian beam, one could draw two conclusions: firstly, sign inconsistence within different frequency components κ(k;ξ) will broaden the Gaussian beam [Figs. 2(c) and 3(f)], which can be seen more clearly in Fig. 1(b), illustrating the dependence of the integral form-factor on the linear chirp. It is worth emphasizing here that there exist two critical values of C: i.e., C₁ ≈ –1.6, below which χ always equals 0.2821 (the input integral form-factor) irrespective of ξ; C₂ ≈ 1.5, above which the Gaussian beam becomes diffraction-free (i.e., χ = 0.1995) after passing the shift point (i.e., ξ ≈ 33 for C = 15). The shift point increases with C. Secondly, once κ (k;ξ) changes its sign, the beam will be “reflected”. However, one needs to bear in mind that the reflection-like behavior presented here is different from the one addressed in Ref. [21]: at the reflection point, the beam in our model is broadened. On the contrary, the beam is sharply compressed as described in Ref. [21].

3.2. Case 2: d(ξ) = cos(Ω ξ)

Under this condition, we also take C = 0 for instance and find that the corresponding trajectory reduces to κ (k;ξ) = k – sin (Ω ξ)/2Ω, which means that the sign of central frequency component (i.e., k = 0) changes at ξ = nξ₀ (where ξ₀ = π/Ω). By the dint of the method used in case 1, the approximate expression for Eq. (11) is given as follows:

It is worth emphasizing here that in the odd half-period ξ ∈ (2mξ₀,2mξ₀ +ξ₀), κ < 0 and ϕ (2mξ₀ → 2mξ₀ +ξ₀) = 1/Ω² – kξ₀; while in the even half-period ξ ∈ (2mξ₀ + ξ₀,2mξ₀ +2ξ₀), κ > 0 and ϕ (2mξ₀ + ξ₀ → 2mξ₀ + 2ξ₀) = 1/Ω² + kξ₀. Thus we can obtain the analytical solution to Eq. (10) in real space as follows:

where A₁ = exp [ –i (4m + 1–cos (Ωξ))/4Ω²], and A₂ = exp [–i (4m + 3 + cos (Ω ξ))/4Ω²]. We can deduce the trajectory of the beam in real space,

From Eq. (19), one can note that the propagation trajectory will exhibit a zigzag pattern, with the turning points locating at ξ = nξ₀, which are also the changing points of the sign of κ. The emerging point of the beam is confined within a limited domain in η spaces (i.e., η ∈ (η₀ – ξ₀/2,η₀)) [Fig. 4(a)]. Again, a good agreement between the above theoretical predictions and numerical results can be illustrated in Figs. 4(a) and 4(b).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Same as Fig. 2, but for d(ξ) = cos(κ ξ) with (a)–(c)C = 0, (d)–(f)C = 5, and (g)–(i)C = –15. In all these panels, κ = 0.05.

For C ∈ (–k_m – 1/(2Ω),k_m +1/(2Ω)), the beam propagates along the pieces of end-to-end straight lines (solid line), which is same as the C = 0 situation and is in full accordance with the theoretical results plotted by dotted line [Fig. 4(d)]. Those turning points (i.e., ξ = 10.5,52.4,136.1,178) are the zeros of the equation C – sin (Ω ξ) / (2Ω) = 0 (i.e., ξ = arcsin(2ΩC/Ω). And the χ (ξ) curves for this two settings share the similar step-like tendency [Figs. 4(c) and 4(f)]. Due to the unsymmetry of κ (ξ) (when the involved chirp is positive (or negative), the corresponding κ (ξ) > 0 (or κ (ξ) < 0) part will account for a larger proportion), the beam gradually deviates from the η = 0 axis [Figs. 4(d) and 4(h)], along which the equivalent refractive index is the largest. It is worth emphasizing is that the diffraction-free beam can also be realized for | C | < 1/(2Ω) + k_m, which ensures that k + C – sin (Ωξ)/(2Ω) ≷ 0 holds regardless of both k and ξ [Figs. 4(g)–4(i)].

3.3. Case 3: d(ξ) = ξ –10

In Fig. 5, we give the comparisons between the theoretical results and the numerical simulation of the evolution behavior in κ and η spaces to study the effect of involved chirp on the dynamics characteristics of a Gaussian beam. According Eq. (9), the trajectory of all the frequency components is

For C ∈(–k_m, + ∞), there is only one zero [Fig. 5(b)], i.e., , so the beam will experience one reflection, and the corresponding trajectory in real space is

which is similar to the d(ξ) = 1 case [comparing Figs. 3(d) and 5(a)]. However, the broadening effect induced by the d(ξ) = ξ –10 modulation is much weaker than that induced by the d(ξ) = 1 modulation (i.e., χ changes from 0.2821 to 0.2812 in the former case, while in the latter case, χ changes from 0.2821 to 0.1995) [Fig. 5(c)]. For C ∈ (–22.5 – k_m,–k_m), two zeros appear [Fig. 5(e)], thus the beam will be reflected twice [Fig. 5(d)]. For the interval C ∈ (–∞,–22.5–k_m), κ (k;ξ) is always smaller than zero irrespective of both k and ξ, which means that the diffraction-free beam is obtained again [Figs. 5(g)–5(i)].

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The same as Fig. 2, but for d(ξ) = ξ – 10 with (a)–(c) C = 15, (d)–(f) C = –15, and (g)–(i) C = –30.

If 1 < α ⩽ 2, which means fractional diffraction effect is involved, the diffraction-free beam will no longer exists [Comparing Fig. 5(g) with Figs. 6(c1)–6(c3)]. From Figs. 6(a1)–6(b3), one can note that the zigzag trajectory propagation still exists. However, the increasement of the Lévy index will result in the transverse drift of those turning points. Taking the d(ξ) = cos (Ω ξ) modulation as an example, as shown in Figs. 6(b1)–6(b3), we have found that the first turning point for different Lévy indices: (62.8,–21.4) for α = 1, (62.8,–59.2) for α = 1.3, (62.8,–140) for α = 1.6, and (62.8,–304) for α = 1.9. In the meantime, as α grows, a clearer delocalization of the Gaussian beam can be seen in Fig. 6.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. The 1 < α ⩽ 2 case. (a1)–(a3) d(ξ) = 1 modulation with C = 15; (b1)–(b3) d(ξ) = cos (Ω ξ) modulation with C = 0; (c1)–(c3) d(ξ) = ξ –10 modulation with C = –30.

In summary, with involving three different kinds of longitudinal modulation (i.e., d(ξ) = 1, d(ξ) = cos (Ωξ), d(ξ) = ξ –10), we have proposed a novel scheme for trajectory engineering of a chirped Gaussian beam modeled by the fractional Schrödinger equation with dynamic linear index potential. For the limiting case α=1, by carefully choosing the linear chirp, the beam will exhibit two totally different evolution characteristics. Firstly, as the whole spatial spectrum lies in the lower half-plane (or upper half-plane), resulting from the frequency shift, the diffraction-free beam propagating along a straight line can be realized. Secondly, if there exist one or more zeros in the spatial spectrum, the beam will experience one or more reflections, which leads to the formation of the zigzag trajectory propagation and the emergence of broadening effect. Importantly, a good agreement between the numerical simulations and the theoretical results has been obtained. For 1 < α ⩽ 2, increasement of the Lévy index not only results in the drift of those turning points along the transverse direction, but also leads to the delocalization of the Gaussian beam. In view of all the above achievements, we believe that our proposal can be applied to other types of beams (i.e., airy beam) and can broaden the potential applications in developing new strategies for trajectory engineering.